home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / TRY / util.f < prev   
Encoding:
Text File  |  1994-08-02  |  9.7 KB  |  432 lines
  1. #define MAX_SIZE 1111
  2. #define MAX_THREADS 8
  3.  
  4. *****************************************************************************
  5. *
  6. *    Array generation
  7. *
  8. *****************************************************************************
  9.       subroutine sgen( a, n)
  10.       real a(*)
  11.       integer n
  12.       integer i,ir
  13.  
  14.       ir = mod(1325 * (irand()+1), 65536)
  15.       ratio = 1.0 / 32768.0
  16.     do i = 1, n
  17.         ir = mod(ir*3125,65536)
  18.         a(i) = float( ir ) * ratio - 1.0
  19.     end do
  20.       return
  21.       end
  22.  
  23.       subroutine dgen( a, n)
  24.       double precision a(*)
  25.       integer n
  26.       integer i,ir
  27.  
  28.       ir = mod(1325 * (irand()+1), 65536)
  29.       ratio = 1.0 / 32768.0
  30.     do i = 1, n
  31.         ir = mod(ir*3125,65536)
  32.         a(i) = float( ir ) * ratio - 1.0
  33.     end do
  34.       return
  35.       end
  36.  
  37.       subroutine cgen( a, n)
  38.       complex a(*)
  39.       integer n
  40.       integer i,ir
  41.       real re, im, ratio
  42.  
  43.       ir = mod(1325 * (irand()+1), 65536)
  44.       ratio = 1.0 / 32768.0
  45.     do i = 1, n
  46.         ir = mod(ir*3125,65536)
  47.         re = float( ir ) * ratio - 1.0
  48.         ir = mod(ir*3125,65536)
  49.         im = float( ir ) * ratio - 1.0
  50.         a(i) = cmplx( re, im)
  51.     end do
  52.       return
  53.       end
  54.  
  55.       subroutine zgen( a, n)
  56.       double complex a(*)
  57.       integer ld, n
  58.       integer i,j,ir
  59.       double precision re, im, ratio
  60.  
  61.       ir = mod(1325 * (irand()+1), 65536)
  62.       ratio = 1.0 / 32768.0
  63.     do i = 1, n
  64.         ir = mod(ir*3125,65536)
  65.         re = float( ir ) * ratio - 1.0
  66.         ir = mod(ir*3125,65536)
  67.         im = float( ir ) * ratio - 1.0
  68.         a(i) = cmplx( re, im)
  69.     end do
  70.     
  71.       return
  72.       end
  73.  
  74. *****************************************************************************
  75. *
  76. *    Array Sum
  77. *
  78. *****************************************************************************
  79.       real function ssum( n, x, incx)
  80.       integer n, incx
  81.       real x(*)
  82.     ssum = 0.0d0
  83.     ix = 1
  84.     do i = 1, n
  85.         ssum = ssum + x(ix)
  86.         ix = ix + incx
  87.     enddo
  88.       return
  89.       end
  90.  
  91.       double precision function dsum( n, x, incx)
  92.       integer n, incx
  93.       double precision x(*)
  94.     dsum = 0.0d0
  95.     ix = 1
  96.     do i = 1, n
  97.         dsum = dsum + x(ix)
  98.         ix = ix + incx
  99.     enddo
  100.       return
  101.       end
  102.  
  103. *****************************************************************************
  104. *
  105. *    1D FFT
  106. *
  107. *****************************************************************************
  108.       subroutine sft1d( plusminus, n, a, inc)
  109.       real a(*)
  110.       real w(MAX_SIZE)
  111.       integer n, plusminus, inc
  112.       
  113.       real u, v
  114.       TPI = 6.28318530717959
  115.  
  116.       if ( plusminus .eq. -1) then
  117.     w(1) = ssum( n, a, inc)
  118.     l = (n+1)/2
  119.     do k = 2, l
  120.         u = 0.0
  121.         v = 0.0
  122.         do i = 1, n
  123.         u = u + a(inc*(i-1)+1) * cos((k-1) * (i-1) * tpi / float(n))
  124.         v = v - a(inc*(i-1)+1) * sin((k-1) * (i-1) * tpi / float(n))
  125.         end do
  126.         w(2*k-2) = u
  127.         w(2*k-1) = v
  128.     end do
  129.     if( mod(n,2) .eq. 0) then
  130.         u = 0.0
  131.         do k = 1, l
  132.         u = u + a(inc*(2*k-1-1)+1) - a(inc*(2*k-1)+1)
  133.         end do
  134.         w(n) = u
  135.     end if
  136.     do i = 1, n
  137.         a(inc*(i-1)+1) = w(i)
  138.     end do
  139.       elseif (plusminus .eq. 1) then
  140.     Print *, 'Sorry NOT Implemented'
  141.       else
  142.     print *, 'Error first argument of SFT1D should -1 or +1'
  143.     stop
  144.       endif
  145.  
  146.       return
  147.       end
  148.  
  149.       subroutine dft1d( plusminus, n, a, inc)
  150.       double precision a(*)
  151.       double precision w(MAX_SIZE)
  152.       double precision dsum
  153.       integer n, plusminus, inc
  154.       
  155.       double precision u, v, tpi
  156.       TPI = 2 * 3.14159265358979323846
  157.  
  158.       if ( plusminus .eq. -1) then
  159.     w(1) = dsum( n, a, inc)
  160.     l = (n+1)/2
  161.     do k = 2, l
  162.         u = 0.0
  163.         v = 0.0
  164.         do i = 1, n
  165.         u = u + a(inc*(i-1)+1) * cos((k-1) * (i-1) * tpi / float(n))
  166.         v = v - a(inc*(i-1)+1) * sin((k-1) * (i-1) * tpi / float(n))
  167.         end do
  168.         w(2*k-2) = u
  169.         w(2*k-1) = v
  170.     end do
  171.     if( mod(n,2) .eq. 0) then
  172.         u = 0.0
  173.         do k = 1, l
  174.         u = u + a(inc*(2*k-1-1)+1) - a(inc*(2*k-1)+1)
  175.         end do
  176.         w(n) = u
  177.     end if
  178.     do i = 1, n
  179.         a(inc*(i-1)+1) = w(i)
  180.     end do
  181.       elseif (plusminus .eq. 1) then
  182.     Print *, 'Sorry NOT Implemented'
  183.       else
  184.     print *, 'Error first argument of SFT1D should -1 or +1'
  185.     stop
  186.       endif
  187.  
  188.       return
  189.       end
  190.  
  191.  
  192. *****************************************************************************
  193. *                                        *
  194. *    2D FFT                                    *
  195. *                                        *
  196. *****************************************************************************
  197.       subroutine sft2d( plusminus, n1, n2, w, ldw)
  198.       real w(ldw,n2)
  199.       real save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
  200.       integer n1, n2, lda, ldw, plusminus
  201.       
  202.  
  203.       if (plusminus .eq. -1) then
  204.     call sfft1di( n1, save)
  205. c$doacross local(j), share(cp)
  206.     do j = 1, n2
  207.         call sfft1d( plusminus, n1, w(1,j), 1, save)
  208.     end do
  209.  
  210.     call sfft1di( n2, save)
  211.     call sfft1d( plusminus, n2, w(1,1), ldw, save)
  212.     if( mod(n1,2) .eq. 0) then
  213.         call sfft1d( plusminus, n2, w(n1,1), ldw, save)
  214.     end if
  215.  
  216.     call cfft1di( n2, save)
  217. c$doacross local(i,j,my)
  218.     do i = 2, n1-1, 2
  219.         my = 1
  220. c$        my = mp_my_threadnum()+1
  221.         do j = 1, n2
  222.         cp(2*j-1,my) = w(i+0,j)
  223.         cp(2*j-0,my) = w(i+1,j)
  224.         end do
  225.         call cfft1d( plusminus, n2, cp(1,my), 1, save)
  226.         do j = 1, n2
  227.         w(i+0,j) = cp(2*j-1,my) 
  228.         w(i+1,j) = cp(2*j-0,my)
  229.         end do
  230.     end do
  231.       else
  232.     Print *, 'Sorry NOT Implemented'
  233.       end if
  234.       return
  235.       end
  236.       
  237.  
  238.       subroutine dft2d( plusminus, n1, n2, w, ldw)
  239.       double precision w(ldw,n2)
  240.       double precision save( 2*MAX_SIZE), cp(2*MAX_SIZE,MAX_THREADS)
  241.       integer n1, n2, lda, ldw, plusminus
  242.       
  243.  
  244.       if (plusminus .eq. -1) then
  245.     call dfft1di( n1, save)
  246. c$doacross local(j)
  247.     do j = 1, n2
  248.         call dfft1d( plusminus, n1, w(1,j), 1, save)
  249.     end do
  250.  
  251.     call dfft1di( n2, save)
  252.     call dfft1d( plusminus, n2, w(1,1), ldw, save)
  253.     if( mod(n1,2) .eq. 0) then
  254.         call dfft1d( plusminus, n2, w(n1,1), ldw, save)
  255.     end if
  256.  
  257.     call zfft1di( n2, save)
  258. c$doacross local(i,j,my), share(cp)
  259.     do i = 2, n1-1, 2
  260.         my = 1
  261. c$        my = mp_my_threadnum()+1
  262.         do j = 1, n2
  263.         cp(2*j-1,my) = w(i+0,j)
  264.         cp(2*j-0,my) = w(i+1,j)
  265.         end do
  266.         call zfft1d( plusminus, n2, cp(1,my), 1, save)
  267.         do j = 1, n2
  268.         w(i+0,j) = cp(2*j-1,my) 
  269.         w(i+1,j) = cp(2*j-0,my)
  270.         end do
  271.     end do
  272.       else
  273.     Print *, 'Sorry NOT Implemented'
  274.       end if
  275.       return
  276.       end
  277.       
  278.  
  279. *****************************************************************************
  280. *
  281. *    3D FFT
  282. *
  283. *****************************************************************************
  284.       subroutine sft3d(plusminus,n1,n2,n3,w,ld1,ld2)
  285.       real w(ld1,ld2,n3)
  286.       real save(2*MAX_SIZE)
  287.       integer n1,n2,n3,ld1,ld2,plusminus
  288.  
  289.       if( Plusminus .ne. -1) then
  290.     Print *, 'Sorry NOT Implemented'
  291.     return
  292.       endif
  293. c
  294. c   transform along X  first ...
  295. c
  296.       call sfft1di( n1, save)
  297. c$doacross local(i,j,k)
  298.       do k = 1, n3
  299.        do j = 1, n2
  300.         call sfft1d( plusminus, n1, w(1,j,k), 1, save)
  301.        end do
  302.       end do
  303. c
  304. c   transform along Y then ...
  305. c
  306.       call sfft1di( n2, save)
  307. c$doacross local(i,j,k)
  308.     do k = 1,n3
  309.         call sfft1d( plusminus, n2, w(1,1,k), ld1, save)
  310.     end do
  311.       if (mod(n1,2) .eq. 0) then
  312. c$doacross local(i,j,k)
  313.     do k = 1,n3
  314.         call sfft1d( plusminus, n2, w(n1,1,k), ld1, save)
  315.     end do
  316.       end if
  317.  
  318.       call cfft1di( n2, save)
  319. c$doacross local(i,k)
  320.       do k = 1,n3
  321.     do i = 2, n1-1, 2
  322.         call csfft1d( plusminus, n2, w(i,1,k), 1, ld1, save)
  323.     end do
  324.       end do
  325. c
  326. c   transform along Z finally ...
  327. c
  328.       call sfft1di( n3, save)
  329.     call sfft1d( plusminus, n3, w(1,1,1), ld1*ld2, save)
  330.     if ( mod(n2,2) .eq. 0) then
  331.         call sfft1d( plusminus, n3, w(1,n2,1), ld1*ld2, save)
  332.     end if
  333.       if (mod(n1,2) .eq. 0) then
  334.         call sfft1d( plusminus, n3, w(n1,1,1), ld1*ld2, save)
  335.     if ( mod(n2,2) .eq. 0) then
  336.         call sfft1d( plusminus, n3, w(n1,n2,1), ld1*ld2, save)
  337.     end if
  338.       end if
  339.  
  340.       call cfft1di( n3, save)
  341. c$doacross local(i,j,k)
  342.     do j = 2, n2-1, 2
  343.         call csfft1d( plusminus, n3, w(1,j,1), ld1, ld1*ld2, save)
  344.         if (mod(n1,2) .eq. 0) then
  345.         call csfft1d( plusminus, n3, w(n1,j,1), ld1, ld1*ld2, save)
  346.         end if
  347.     end do
  348. c$doacross local(i,j)
  349.       do j = 1,n2
  350.     do i = 2, n1-1, 2
  351.         call csfft1d( plusminus, n3, w(i,j,1), 1, ld1*ld2, save)
  352.     end do
  353.       end do
  354.  
  355.       return
  356.       end
  357.  
  358.       subroutine dft3d(plusminus,n1,n2,n3,w,ld1,ld2)
  359.       double precision w(ld1,ld2,n3)
  360.       double precision save(2*MAX_SIZE)
  361.       integer n1,n2,n3,ld1,ld2,plusminus
  362.  
  363.       if( Plusminus .ne. -1) then
  364.     Print *, 'Sorry NOT Implemented'
  365.     return
  366.       endif
  367. c
  368. c   transform along X  first ...
  369. c
  370.       call dfft1di( n1, save)
  371. c$doacross local(i,j,k)
  372.       do k = 1, n3
  373.        do j = 1, n2
  374.         call dfft1d( plusminus, n1, w(1,j,k), 1, save)
  375.        end do
  376.       end do
  377. c
  378. c   transform along Y then ...
  379. c
  380.       call dfft1di( n2, save)
  381. c$doacross local(i,j,k)
  382.     do k = 1,n3
  383.         call dfft1d( plusminus, n2, w(1,1,k), ld1, save)
  384.     end do
  385.       if (mod(n1,2) .eq. 0) then
  386. c$doacross local(i,j,k)
  387.     do k = 1,n3
  388.         call dfft1d( plusminus, n2, w(n1,1,k), ld1, save)
  389.     end do
  390.       end if
  391.  
  392.       call zfft1di( n2, save)
  393. c$doacross local(i,k)
  394.       do k = 1,n3
  395.     do i = 2, n1-1, 2
  396.         call zdfft1d( plusminus, n2, w(i,1,k), 1, ld1, save)
  397.     end do
  398.       end do
  399. c
  400. c   transform along Z finally ...
  401. c
  402.       call dfft1di( n3, save)
  403.     call dfft1d( plusminus, n3, w(1,1,1), ld1*ld2, save)
  404.     if ( mod(n2,2) .eq. 0) then
  405.         call dfft1d( plusminus, n3, w(1,n2,1), ld1*ld2, save)
  406.     end if
  407.       if (mod(n1,2) .eq. 0) then
  408.         call dfft1d( plusminus, n3, w(n1,1,1), ld1*ld2, save)
  409.     if ( mod(n2,2) .eq. 0) then
  410.         call dfft1d( plusminus, n3, w(n1,n2,1), ld1*ld2, save)
  411.     end if
  412.       end if
  413.  
  414.       call zfft1di( n3, save)
  415. c$doacross local(i,j,k)
  416.     do j = 2, n2-1, 2
  417.         call zdfft1d( plusminus, n3, w(1,j,1), ld1, ld1*ld2, save)
  418.         if (mod(n1,2) .eq. 0) then
  419.         call zdfft1d( plusminus, n3, w(n1,j,1), ld1, ld1*ld2, save)
  420.         end if
  421.     end do
  422. c$doacross local(i,j)
  423.       do j = 1,n2
  424.     do i = 2, n1-1, 2
  425.         call zdfft1d( plusminus, n3, w(i,j,1), 1, ld1*ld2, save)
  426.     end do
  427.       end do
  428.  
  429.       return
  430.       end
  431.  
  432.